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1 Introduction 

The discovery of planets orbiting in binary star systems represents an exciting 
new field of astrophysics. The stability of planetary orbits in binary systems 
can only be addressed analytically in special cases, so most researchers have 
studied stability using long-term N-body integrations of test particles, exam- 
ining binary systems with a range of masses and orbits (e.g. Wiegert and 
Holman 1997, Holman and Wiegert 1999, Haghighipour 2006). This has led 
to a good understanding of the likely regions of stability and instability in 
binary systems. Integrators can also been used to study the more complex 
problem of several finite-mass planets orbiting in a binary system, where in- 
teractions between the planets are significant. However, at the time of writing, 
this problem has been explored in less detail than the test-particle case, and 
we still lack a general theory for the stability of these systems. 

N-body integrators have found a second application modelling the forma- 
tion of planetary systems around binary stars (e.g. Quintana et al. 2006). 
Typically, these studies have paralleled those of planet formation around sin- 
gle stars, examining a particular stage of growth such as the formation of 
planetesimals, oligarchic growth, or late-stage accretion of terrestrial planets. 
The results of these studies are discussed extensively in the chapter in this 
volume by Quintana et al. 

Most conventional integrator algorithms can be applied to binary star 
systems with little or no modification. Runge-Kutta, Bulirsch-Stoer and Ev- 
erhart's RADAU integrators fall into this category for example (Press et al. 
1992, Stoer and Bulirsch 1980, Everhart 1985). These algorithms contain no 
built-in information about the system of differential equations they are solv- 
ing, so they can be applied to binary systems and single-star systems equally 
well. 

Over the last decade and a half, symplectic integrator algorithms have 
become increasingly popular and are widely used to study the dynamics of 
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planetary and satellite orbits. These algorithms have two advantages over con- 
ventional integrators. First, symplectic integrators typically have good long- 
term energy conversation properties. While energy is not conserved in most 
problems, the energy error typically makes high frequency oscillations about 
zero, while exhibiting no long-term trend beyond that generated by computer 
round-off error. Secondly, in problems involving a dominant, primary mass, 
such as the Sun in the Solar System, the motion of other objects about the 
central body can be "built in" . A relatively small amount of computation is 
required to calculate the accelerations due to the central body, and a large 
stepsize can be used since this only needs to be small enough to resolve the 
perturbations between the smaller bodies. These advantages mean that sym- 
plectic integrators have become the tool of choice for many researchers at 
present, and they will be the focus of this chapter. 

The rest of the chapter is organized as follows. Section 2 contains a review 
of the original mixed- variable symplectic mapping developed by Wisdom and 
Holman (1991). Section 3 shows how this algorithm has been modified for use 
specifically in binary-star systems. Section 4 shows how symplectic integrators 
can be improved by developing symplectic correctors, and describes a new 
corrector for binary-star algorithms. Section 5 discusses problems that can 
arise when planets come close to one or both binary stars, and what might be 
done to overcome these problems. Finally Section 6 contains a short summary. 

2 Mixed- Variable Symplectic Integrators 

The most widely used symplectic integrators applied to planetary systems are 
"mixed-variable symplectic" (MVS) mappings, so called because they sepa- 
rate a problem into two parts, each of which is solved using a different set 
of variables (typically Cartesian coordinates and orbital elements). These al- 
gorithms were first introduced by Wisdom and Holman (1991) and described 
independently by Kinoshita et al. (1991) 

To understand how these integrators work, it is easiest to start by consid- 
ering Hamilton's equations for a system of N bodies: 



where x and p are the coordinates and momenta of the bodies respectively, 
and H is the Hamiltonian of the system. Using Hamilton's equations, the 
evolution of any quantity q can be expressed as 



dxi 
~~dt 

dpi 
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(1) 
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where {,} are Poisson brackets, and F is an operator that depends on the 
Hamiltonian. The evolution of q can be found by solving (f5]) , which gives 

q(r) = e TF q(Q) = (l + tF + ^ + • • q(0) (3) 

MVS integrators divide the Hamiltonian into several parts each of which 
can be solved efficiently in the absence of the others. Most algorithms divide 
H into parts that can be solved analytically although this is not strictly nec- 
essary. If we separate the Hamiltonian so that H = Ha + Hb , with operators 
A and B corresponding to the Hamiltonians H a and Hb, then 

q(r) = e^ A +%(0) (4) 

where 

= 1 + r{A + B) + rHJ ± AB + BA ± B^ 1+ ^ 

In general, the operators A and B will not commute so that AB ^ BA. 

A symplectic integrator is generated by concatenating several terms of the 
form exp(aferA) and exp(&fcT.B), where the a& and bk are constant coefficients. 
The goal is to make the resulting expression equal to ([5]) up to some order in 
the stepsize r. This is most easily accomplished using the Baker-Campbell- 
Hausdorff (BCH) formula which expresses the product of two exponential 
operators as a single new exponential operator: 

e A e B =e X p!.A + B+^[A,B} + ^[A,A,B} + ^[B,B,A}+--\ (6) 

where [A, B] = AB - BA and [A, B, C] = [A, [B, C]] etc. (Yoshida 1990). 

The most commonly used MVS algorithm is the second-order leapfrog 
integrator: 

exp (Ia) exp(rB) exp (Ia) = exp Sr(A + B) + ^[B, B, A] 

- £[A,A,B]+0(t*)\ (7) 



This differs from the true evolution ([5]) by terms proportional to r 3 and higher. 
Over the course of a long integration, the number of steps will be inversely 
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proportional to r. The total error will be 0(t 2 ), meaning that leapfrog is a 
2nd-order integrator. 

Each timestep using the leapfrog algorithm consists of advancing the sys- 
tem corresponding to Ha for a time r/2, then advancing Hb for r, and finally 
advancing Ha for r/2. The integrator consists of 3 substeps. However, the first 
and last of these both involve A, so the last substep of one timestep can be 
combined with the first substep of the following step. Over the course of a 
long integration, each timestep is effectively composed of only 2 substeps. 

Wisdom and Holman (1991) split the Hamiltonian into a part i?Kc P con- 
taining terms corresponding to the Keplerian motion of each planet about the 
central star, and a second part -ffint containing direct and indirect perturba- 
tion terms due to interactions between the planets. This is accomplished using 
Jacobi coordinates, where the position of the innermost planet is measured 
with respect to the central star, and the positions of the remaining planets are 
measured with respect to the centre of mass of the central star and planets 
with lower indices. Evolution under the Keplerian part of the Hamiltonian 
can be calculated efficiently using Gauss's / and g functions (Danby 1988). If 
Jacobi coordinates are used, H\ nt is a function of the coordinates only, so this 
part of the problem can also be solved analytically. Jacobi coordinates are 
also the natural choice for hierarchical systems like the planets in the Solar 
System, and generally lead to smaller errors than other canonical coordinate 
systems, such as barycentric coordinates, for a given step size. Barycentric 
coordinates are an especially poor choice for the Solar System since the inner 
planets are less massive than the outer ones. As a result, the guiding centre for 
the motion of the inner planets is much closer to the Sun than the barycentre 
of the Solar System. 

When applied to planetary systems, H\ nt ~ ei?Ke P in the Wisdom-Holman 
mapping, where e is the planetary to stellar mass ratio, which is typically 
small. The error over a long integration is therfore 0(er 2 ), and the small 
value of e ensures that MVS leapfrog performs well even though it is only a 
second-order integrator. 

3 Binary-Star Algorithms 

The mixed-variable symplectic (MVS) integrators described in the previous 
section can be applied to any hierarchical system of bodies. This means they 
can be used to calculate the orbital evolution of planets in a binary star sys- 
tem provided that the radial ordering of the planets doesn't change. The MVS 
algorithm can also be adapted to systems that contain multiple hierarchies, 
such as two binary systems in orbit about each other, by using a generalized 
version of Jacobi coordinates (Beust 2003). However, whenever planets come 
close to one another, the condition Hi nt <C Hk cp is violated and MVS inte- 
grators will perform poorly. In situations where planets have eccentric orbits 
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that cross those of their neighbours, or where close encounters are possible, 
the MVS algorithms must be modified. 

Duncan et al. (1998) and Chambers (1999) have described two ways to 
do this. Here we will describe the latter approach since this leads directly to 
a new class of symplectic integrators that can be applied to orbits in binary 
systems. For reasons that will become clear, the new method requires a new 
set of coordinates. Ideally, these should include three spatial coordinates for 
the centre of mass of the system (as Jacobi coordinates do), and treat all the 
planets cquivalently, that is, make no asumptions about their radial ordering 
(in contrast to Jacobi coordinates). 

Canonical heliocentric coordinates (also called democratic heliocentric co- 
ordinates) meet both of these requirements (Duncan et al. 1998). Here, the 
position of each planet is measured with respect to that of the central star, 
and the stellar coordinates are replaced with those of the centre of mass of 
the system: 

m x + E^Li rtij-x-j 
X = 

TOtot 

Xj = Xj - x (8) 

where subscript refers to the star, rm is the mass of planet i, and m to t is 
the total mass of the system. 

The canonically conjugate momenta (which correspond to barycentric ve- 
locities) are: 



JV 

Po = Po + X] 



m. 



Pi = Pi-— *-(!*> + 5>i) ( 9 ) 

Using these coordinates, the Hamiltonian for a system of TV planets orbit- 
ing a single star can be split into three parts: 



H — i?Kcp + -fflnt + ffjump (10) 



where 



JV 



I p i Gm m, 



1=1 

JV 



2rrii Ri 



\ \^ Gmiirij 



H 



Jump 



==(x> <"» 



,i=l 
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where Ri = |X, | and = |Xj — Xj|. Note that we have dropped a term 
Pq /2m tot which simply acts to move the centre of mass at a constant velocity. 

Several second-order symplectic integrators can be constructed using 
canonical heliocentric coordinates, for example: 



where /, J and K are operators associated with Hi nt , Hj ump and H^cp respec- 
tively. Other second-order algorithms are similar except that the operators are 
permuted, making sure that the arrangement is symmetrical in each case. 

Advancing the system under H\ nt is straightforward since this part of the 
Hamiltonian is a function of the coordinates only. As a result, the positions 
of all planets remain constant while the velocities change due to perturba- 
tions from the other planets. Advancing under i?j um p is trivial since this is 
a function of the momenta only. In this case, each planet's velocity remains 
constant but its spatial coordinates jump by a small amount. This jump is the 
same for all planets, so it becomes relatively more important for objects close 
to the central star, a point we will return to in Section 5. Advancing i?Kep is 
best done using Gauss's / and g functions as before, noting that and X.; 
are canonically conjugate. 

The integrator described by (fT!2"|) performs quite well for planetary systems 
in which the planets do not undergo close encounters. Typical energy errors 
are intermediate between those using Jacobi coordinates and barycentric co- 
ordinates for hierarchical systems, although Jacobi coordinates lose their ad- 
vantage if the planets have crossing orbits. Despite consisting of 5 substeps 
rather than 3, the integrator described above involves slightly less computa- 
tional effort than the MVS leapfrog algorithm since advancing under Hj ump is 
trivial, while Hi nt only contains direct terms whereas indirect terms are also 
present when using Jacobi coordinates. 

As with the MVS leapfrog integrator, problems arise when a pair of plan- 
ets has a close encounter because -Hint can become comparable in size to 
Hkc-p- Chambers (1999) showed that this difficulty can be overcome using a 
hybrid algorithm. Here, each term in Hi nt is split between Hi Dt and Hk cp so 
that the former always remains much smaller than the latter. Under this new 
arrangement, the Hamiltonian is divided as follows: 




(12) 





(13) 
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where r is a partition function. 

A second-order hybrid integrator has the same form as (|12|) . that is: 

exp \y) exp (t) exp ^ tV> exp (t ) exp (t) ^ 

where L and S are operators associated with H^ aTge and -HsmaU respectively. 

The partition function is chosen so that r(R) = 1 when R is large and 
r(R) — > as i? — > 0. With this choice of J 1 , it is always the case that 
H Large ^ -^Smaii, and the resulting integrator remains accurate during close 
encounters between planets. However, -f/Large is no longer analytically soluble 
during a close encounter since it now includes a three-body problem. These 
three-body terms can be calculated using a conventional TV-body algorithm 
such as Bulirsch-Stoer. Provided this is done to an accuracy level close to 
machine precision, the user shouldn't be able to tell in practice whether the 
solution was derived analytically or numerically. 

Using Bulirsch-Stoer (or any other strategy) to evolve the system through 
a close encounter will slow down an integration. However, only the pair of 
planets involved in the encounter need to be integrated in this way. All the 
other planets are advanced analytically under i/Largc using Gauss's functions. 
When there are no encounters, the algorithm becomes identical to (fT2|) and 
there is no loss of speed. 

The algorithm given by (|14() works well for most systems of planets orbiting 
a single star, and has also been used extensively to study the formation of 
planets from a disk of smaller bodies, since these bodies undergo many close 
encounters. However, the algorithm performs poorly when applied to planets 
orbiting in binary systems. The problem arises because Hj nmp now contains a 
momentum contribution from one of the binary stars. The stellar momentum 
is typically large and this leads to large changes in position for all the planets 
via -ffjump- It is no longer true that Hj ump <C i?Kop, and the error per step 
becomes large unless an unacceptably small timestep is chosen. 

As Chambers et al. (2002) have shown, the solution to this problem is to 
devise new coordinate systems for binary systems such that all large terms 
can be incorporated into a single part of the Hamiltonian. Stable planetary 
orbits typically fall into one of two classes: (i) those that are tightly bound 
to one member of a binary, or (ii) those that orbit both stars at a distance 
that is considerably larger than the semi-major axis of the binary orbit. Each 
configuration will require a different set of coordinates and we will consider 
the two cases separately. 



3.1 Wide Binary Case 



The Hamiltonian for a system containing N planets orbiting one member of 
a binary star is 
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'i=l i=l ]>i J 

where the planets orbit star A while star B is a distant companion. 

Making use of the hierarchical arrangement of the binary system, we define 
a new set of coordinates, called wide-binary coordinates, as follows: 



V m tot J 



X, = Xj — x^4, 

( m A x A + g mpcj \ 

X B = x B - , (16) 

y m A + 22 j m i J 

where m to t — mA + m B + 2~2j frij is the total mass of the system, and each 
of the summations run from 1 to N . Using these coordinates, the position of 
each planet is measured with respect to star A, while the position of star B 
is measured with respect to the center of mass of all the other objects. 
The conjugate momenta P are: 




Pb = Pb - m B [ J - , (17 

V m tot / 

where the summations run from 1 to N. 

In terms of the new coordinates, the Hamiltonian can be written as 

H = Hjicp + Hi n t + -ffjump, (18) 



where 



Kcp ~ I 2/z bin ~ R B J ^2m, Ri 

N 



Hint = -2^2^ R . + Gm BTfl A [ — 

i=l j>i % 3 ^ 

N ( 1 

+ Gm B X! m M "TJ 



i=i 



Rb |Xb — Xj + SI 
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where /ibm = {jnA rn i) m B /^tot is the reduced mass of the binary system 
(including the mass of the planets) , and 

« - ^ i=i ™; 1 (20) 



in 



The terms in -ffKcp consist of those due to the Keplerian motion of the 
binary (adding the masses of the planets to star A) and those due to the 
Keplerian motion of the planets about star A. The terms in Hi nt represent the 
interactions between planets and also the tidal perturbations on the planets 
due to star B. Finally, Hjump contains indirect perturbation terms. 

In the absence of close encounters, Hi nt <C i?Kcp and -ffj ump -C -ffjccp- Each 
part of the Hamiltonian can be advanced efficiently using analytic solutions. 
For example, the x component of the acceleration of planet k under -f/r n t is 
given by 

dV x ,k 1 dHi nt 
dt TOfe dXk 

( GmArriB \ X B + S x 



\mA + 'E,i m i/ |X B +S 



Gm B \ X B - Xi + S x 



mA + EimJ^i t |X B -X. i + S| 3 

|x B -x fc + s| 3 ^ Rf k 

where V is the velocity. Note that the acceleration on planet k does not involve 
any terms proportional to l/m^, so test particles can be integrated in exactly 
the same way as massive planets. 

The x component of the acceleration on star B is given by 



dV x , B _ 
— — = Gm A 
at 



X B X B + S x 



R% |X B + S|3 



N 



m i 



X B (X B — Xi 



R% IXB-X. + SI 3 



(22) 

Close encounters between planets can be dealt with in the same way as 
for systems with a single star by partitioning the planet interaction terms 
between i/xop and Hi nt as in (|13p . 

One step of the new wide-binary algorithm consists of 5 substeps: 

• Advance H\ nt for r/2, where r is the timestep. 

• Advance -Hj ump for r/2. 

• Advance H^ ep for r. 
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• Advance -ffj um p for t/2. 

• Advance Hi nt for t/2. 

The first and last substeps can be combined into a single substep except at 
the beginning of the integration or whenever output is required. 

The wide-binary integrator is a second-order algorithm since the 3 pieces 
of the Hamiltonian are applied in a symmetric order (see Yoshida 1990). Each 
timestep has an error 0(er 3 ), where e is the ratio of the planetary mass to 
the stellar mass so that e<l. 

Figure 1 compares the accuracy of the wide-binary algorithm with the 
hybrid symplectic integrator of (|14p when integrating the four giant planets 
of the Solar System in the presence of a binary companion. The giant planets 
orbit the Sun, while a second solar-mass star orbits the combined system 
moving on an orbit with semi-major axis a = 160 AU, eccentricity e = 0.25 
and inclination of 0. The figure shows the energy error as a function of time for 
a 100,000-year integration using a stepsize of 50 days. The upper panel shows 
the performance of the hybrid integrator, while the lower panel shows the 
wide-binary algorithm. The hybrid algorithm performs poorly since it treats 
the binary companion as the equivalent of an additional planet, so that Hj nmp 
is no longer small compared to H^ cp . The wide-binary algorithm, which treats 
the binary companion as a special body, has an energy error about three orders 
of magnitude lower as a result. 



3.2 Close Binary Case 

The Hamiltonian for a system of N planets orbiting a close binary has the 
same form as (|15|) except that now it is understood that the planets orbit 
both members of the binary. The hierarchical nature of the system suggests 
we switch to the following close-binary coordinates: 

m A Tt A + J2j TO i x j + ™bXb 
= , 

TOtot 

X* = Xj - {va^a + vbxb) , 

X B = x B - x.a, (23) 

where mtot is the total mass of all the bodies, the summations run from 1 to 
N, and 

va = m A /m hin 

v B = m B /m hin (24) 

where Wbm = fn A + tub is the mass of the binary, 

Using these coordinates, the position of each planet is measured with re- 
spect to the center of mass of the two stars, while X# is the relative coordi- 
nates of the binary itself. 
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The conjugate momenta P are 



N 

Pa = PA + PB + ^2pj, 

/PA + PB+EjPi 

P, = p 4 - mA J - 

V "Hot 

Pb = Pb (Pa + Pb) , 

where summation indices run from 1 to N . 
The new Hamiltonian is 

H = H Kcp + Hi nt + -ffj um p, 



(25) 



(26) 



where 



#Kcp — 



H 



Jump 



P B Gm hin fi Vv 



2/Xbin Rb 

N 

H Int = Gm hin ^2 771 



x -> x Grriimj 



E 

i=i l 



2rrii Ri 



vb 



Ri IX. + z/bXbI \Xi-iy A Xi 



(27) 



where /Xbin = tuattib / {mA + "is) is the reduced mass of the binary. 

Using the close-binary coordinates, terms in H^cp correspond to the Ke- 
plerian motion of the two binary stars about their common centre of mass, 
and of the planets about this centre of mass. In addition, H\ nt contains terms 
due to interactions between the planets, and perturbations on the planetary 
orbits caused by higher order moments of the binary potential. As before, 
ffj ump contains indirect correction terms. 

In the absence of close encounters, Hi nt and Hj ump are small compared to 
f/Kep, and each part of the Hamiltonian can be advanced analytically. The x 
component of the acceleration on planet fc, when advancing H\ nt , is given by 



dVx,k _ Gm hin X k _ ^ Gruj ^ _ ^ 



dt 



- G 



ik 

vbX b ) 



(X k - vaXb) 
X fe - ^aX b | 3 



(28) 



while the corresponding acceleration on star B is given by 
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dV, 



x,B 



dt 



G E 



rrii 



{X t - v A X B ) (X, + ubXb) 



(29) 



where Pb — ^binVe 

Close encounters between planets can be included in the same way as 
before by dividing the planet interaction terms between H^cp and i?i n t- 

One could devise a second-order scheme using close-binary coordinates 
that is analogous to the second-order wide-binary integrator described above. 
This scheme would contain 5 terms arranged symmetrically. However, this 
scheme would have to use a small stepsize in order to accurately integrate the 
orbit of the binary star. It is more efficient to assign a separate small stepsize 
T/-/Vbin to the binary star, and choose a larger global sizestep r to integrate 
the planets. (This is analogous to the individual-timestep procedure described 
by Saha and Tremaine 1994.) We can do this by splitting i?Kop into a part 
-^BKcp that involves terms in and Pb and a part ffpKcp that does not. 
In a similar manner, we split H\ nt into 2 new parts if B i n t and ffpint, where 



-ffBKep 
ffpKcp 



P B GmbinMbi) 



2^bi] 

N 



Rb 



= E 

i=l 



Pf _ Gm hin mj 
2m,i Ri 

N 

if Bint = Gmbin m 



1=1 



1 

Ri 



VB 



X, + Vb^I 



IX, - i/^Xi 



N 



Hi 



Pint 



EE 

i= 1 j>i 



GmiUij 
Ri 



(30) 



An efficient second-order close-binary scheme has the following form: 

• Advance i?pi n t for r/2, where t is the timestep. 

• Repeat the following iVbi n times: 

Advance ff B int for r/(2Abin)- 
Advance ff B Ko P for r/(2Abin)- 

• Advance ffj U m P for r/2. 

• Advance ffpKop for r. 

• Advance ffj U m P for r/2. 

• Repeat the following iVbi n times: 

Advance ff B Ko P for r/(2Abin)- 
Advance ff B int for T/(2A^ b in)- 

• Advance ffpint for r/2. 

Chambers et al. (2002) suggested making iVbi n smaller than the global 
timestep by a factor equal to the ratio of the binary orbital period to the 
period of the innermost planet. One could also use individual timesteps for 
the binary companion and the planets in the wide-binary algorithm described 
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earlier. However, the amount of computer time saved would be modest since 
most of the effort is required to calculate the direct perturbations between 
the planets and this would not change using individual timesteps. 



4 Symplectic Correctors 

Wisdom et al. (1996) showed that the performance of the MVS mapping can 
be improved at little extra cost by applying "symplectic correctors" . If exp(Af) 
represents a single step of an integrator, given by |7]) for example, the addition 
of a corrector modifies the step to become 

e c e M e- c (31) 

where C is chosen in order to remove the leading order error terms. Over 
the course of multiple timesteps, the corrector and inverse corrector terms 
cancel out, so these only need to be applied at the beginning and end of an 
integration and when output is required. Thus, the addition of a corrector 
involves little extra computational expense over a long integration. 

Following Yoshida (1990), one uncorrected step of the MVS leapfrog inte- 
grator can be expressed as 

m { tI \ i ts\ ( tI 

e = exp — exp(rii) exp 



= exp |r(if + I) + — [K, K,I]-— [K, K, K, K, I] 

+ 0{er 7 ) + 0(e 2 r 3 )} (32) 

where K and I are operators associated with H^cp and Hint respectively, 
e ~ I/K is the planet to star mass ratio, and we have explicitly listed only 
communitator terms where I occurs once. 

We wish to devise a corrector that eliminates terms in [K, K, I] and 
[K, K, K, K, I] . In this way we can reduce the error per step from the usual 
0(er 3 ) to 0(e 2 r 3 ). 

Using the identity 

e c e M e- c = exp j M + [C, M] + ^ [C, C, M] + ^ [C, C, C, M] -\ j (33) 

we find that a corrector of the form 

e c = exp { ar 2 [K, I] + br 4 [K, K, K, I] + 0(er 6 ) + 0(e 2 T 3 )} (34) 



where a and b are constants, results in a corrected step 
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e c e M e- c = exp 1{K + I)t + - a ) t 3 [K. K. I 



^ + bj r 5 [K, K, K, K, I] + 0(er 7 ) + 0(e 2 r 3 ) | 



(35) 



so that the [K, K, I] and [K, K, K, K, I] terms can be eliminated if a = 1/12 
and b= -1/720. 

To be useful in practice, we need to be able to express the corrector as the 
product of terms involving exp(i^) and exp(J) separately. Following Wisdom 
et al. (1996), correctors of the form (fM)) can be developed by noting that 

Y(i,k) =e kK e u e~ 2kK e- u e kK 

( ik 3 1 

= exp i 2ik[K, I] + — [K,K, K,I] + ■ ■ \ (36) 

where we retain only commutators in which / appears once. 

Combining two such expressions, we obtain a suitable series of operators 
for the corrector (f3~4"|) : 

Y{i x , ki) ■ Y(i 2 , k 2 ) = cxp{2(z 1 fc 1 + i 2 k 2 )[K, 1} 

+ ^{iikl+i 2 k 3 )[K,K,K,I} + --\ (37) 



where to match (jMJ), we require that 



ilki + i 2 k 2 = 2^ 

i^l + i^ 3 (38) 



There are many possible solutions to (|38|) , for example 

\A0 



h 

?2 



72 
3^10 
10 

24 



h - ^ (39) 
5 

which generates the symplectic corrector included in the Mercury integrator 
package (Chambers 1999). Wisdom et al. (1996) provide equivalent and higher 
order correctors for the alternative second-order mapping 

e M = exp (if) exp(rJ) exp (if) (40) 
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Symplcctic correctors can also be devised for the binary-star integrators 
described in Section 3. The problem is different in that the Hamiltonian con- 
sists of three parts rather than two, but this only complicates things slightly 
provided we want a corrector that eliminates only terms O(e) . We start with 
an expression for one step of the wide or close-binary algorithms: 

e M = exp j exp (^pj cxp(rK) exp (^pj CX P ( y 

{3 3 
r (I + J + K) + T - [K, K, I] + T - [K, K, J] 

5 5 ^ 

- ^ {K, K, K, K,I}-^ {K, K, K, K, J] + 0(er 7 ) + O(eV) j (41) 

where I, J and K are operators associated with Hi nt , Hj ump and Hk cp for the 
wide or close-binary integrator, and we list only commutators that contain I 
or J once. 

We wish to eliminate the leading-order error terms involving [K, K, I] and 
[K, K, J] , as well as [K, K, K, K, I] and [K, K, K, K, J] . This suggests we look 
for a corrector of the form 

e c = exp { gr 2 [K, I] + hr 2 [K, J] + pr 4 [K, K, K, I] 

+ qT 4 [K,K,K,J]+0(eT 6 ) + 0(e 2 T 3 )} (42) 

where g, h, p and q are constants, which gives a corrected step of the form 



exp | (7 + J + K)t + (j2 " g ) T [K > K > I] + { U ~ k ) T [K > K > J] 

~ (j^+P)i*[K,K,K,K,I]- [-L +q ys[K,K,K,K,J] 

+ 0(er 7 ) + 0(e 2 r 3 )} (43) 

so that terms in [K, K, I] , [K, K, J] , [K, K, K, K, I] and [K, K, K, K, J] can 
all be eliminated if we choose g = h = 1/12 and p = q = —1/720. 

Following the same procedure used to obtain the MVS corrector, we note 
that 

e c e A/2 e B e A ' 2 e- c = exp j A + B + [C, A] + [C, B] + C, A] + 1[C, C, B] 



+ ±[C,C,C,A] + ±[C,C,C,B] + ---} 



(44) 



where we give only commutators that contain A or B once. 
Hence 

Z(i,j,k) = e kK ei J l 2 e u e3 J l 2 e- 2kK e-i J l 2 e- u e-i J / 2 e kK 
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= exp ^2ik[K, I] + 2jk[K, J] + — [K, K, K, I] 

+ J Jf-[K,K,K,J} + ...} (45) 

retaining only commutators that contain I or J once. 
Combining two such expressions gives 

Z(ii,ji,ki) ■ Z(i 2 ,j 2 ,k 2 ) = cxp{2(iifci +i 2 k 2 )[K,I] 
+2(jifci +j 2 k 2 )[K, J] + ^{hkl+i 2 kl)[K,K, K,I] 

+ \(jik\ + j 2 k 3 2 )[K, K,K,J} + -- | (46) 

To provide a suitable corrector, we need to satisfy the following criteria 

iifci + i 2 k 2 = 
jiki+j 2 k 2 = 
Hkl + i 2 kl = 

nkl+nkl^-^ (47) 

so clearly one possible solution is 



*i = 3i 

h = 



1-2 = ]2 



72 

3^10 
10 



24 



fa = ^ (48) 

Figure 2 shows the effect of including this corrector when rerunning the 
integration shown in Figure 1 . The upper panel shows the energy error versus 
time when integrating the four giant planets of the Solar System with a binary 
companion using the corrector. The lower panel shows the case without a 
corrector — the same as in Figure 1. The corrector reduces the energy error 
by roughly three orders of magnitude, comparable to the Jupiter/Sun mass 
ratio. However, because the corrector is only applied at the beginning of the 
integration and prior to each output, this improvement in accuracy is achieved 
at little computational cost. 
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5 Stellar Encounters 

The binary algorithms described in the previous sections are designed to work 
in particular circumstances. In the close-binary algorithm, the planets are 
assumed to orbit the centre of mass of a binary at a distance large compared 
to the binary separation. If a planet comes sufficiently close to the binary 
stars, this assumption will no longer be valid and the algorithm will break 
down. In the wide-binary algorithm, planets are assumed to orbit one member 
of a binary, while receiving small perturbations from the other star. If the 
distance between a planet and its central star ever becomes comparable to the 
distance to the other star, the wide-binary algorithm will also break down. 
The accuracy of each of the algorithms depends on the hierarchy of the system 
being preserved. For this reason the binary integrators are unable to follow 
the trajectories of planets moving on transfer orbits, where the centre of a 
planet's motion switches from one star to another or from one star to both 
stars. 

The wide-binary algorithm suffers from a second limitation in that a planet 
cannot travel too close to its central star either. If this happens, the fixed 
stepsize of the integrator will be too large to properly follow the planet's 
periastron passage and accuracy will be lost. This is a well known limitation 
of symplectic integrators in general, including the MVS mapping (Rauch and 
Holman 1999). Accuracy can be restored by regularizing the motion, so that 
the stepsize effectively depends on a planet's distance from the star (e.g. Preto 
and Tremaine 1999), but regularization becomes highly inefficient for problems 
involving more than a few bodies. 

Levison and Duncan (2000) have suggested an alternative solution which 
involves a new division of the Hamiltonian such that a planet's indirect per- 
turbation terms are added to the Keplerian part of the Hamiltonian whenever 
these terms become large. This is analogous to the procedure described ear- 
lier for maintaining accuracy during a close encounter between two planets. 
We can see how this works by considering the Hamiltonian for the single-star 
case described by (1111) . In Levison and Duncan's scheme, the terms in -ffjump 
are divided between Hj nmp and i?Kep in such a way that the former always 
remains small compared to the latter: 



Levison and Duncan (2000) advocate using a partition function of the form 






(49) 
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JV 

A = 1 - JJ[1 - HRi)] (50) 

i=l 

where X(R) is chosen so that A—*l when any planet approaches the star and 
A — when all the planets are far from the star. As with the hybrid integrator 
described by (flU)) . H^ algc has to be integrated numerically whenever A ^ 0. 
In addition, for this choice of A, Hsmall must also be integrated numerically. 
An obvious shortcoming of this procedure is that an integration will proceed 
slowly whenever any planet passes close to the star since all the planets have 
to be integrated numerically in this case. However, if close periastron passages 
are relatively rare, this shortcoming is not severe. 

The scheme of Levison and Duncan (2000) can be applied to the wide- 
binary algorithm to improve the accuracy whenever a planet passes close to 
the central star. A similar procedure could be developed to cope with planets 
that stray far from the central star in a wide binary, or come close to the stars 
in a close binary. In either case, the planet's orbit is likely to be unstable, 
so this state of affairs will be short-lived, compensating for the fact that 
all objects will have to be integrated numerically during this stage of the 
evolution. 

However, this approach to dealing with stellar encounters suffers from 
two other drawbacks that limit its usefulness. The fact that -ffs ma ii must be 
integrated numerically in addition to -ffLargc can be overcome by choosing a 
partition function X(V) that depends on the planets' velocities rather than 
their positions. Since a planet's velocity will typically become large whenever 
it approaches a star, velocity can be used instead of position to identify a close 
encounter with a star. 

A more serious problem arises when one considers low-mass planets or test 
particles. When advancing the system under Hsmall) the rate of change of the 
x component of velocity of planet k is given by Hamilton's equations: 

AVx.k _ 1 9-ffsmall 

dt mk dXk 

Note that the righthand side of this expression is proportional to 1/rrik, 
so the rate of change of the planet's velocity will become large if the planet's 
mass is small, and will become infinite for massless test particles. This means 
the scheme cannot be used to integrate test particles, and the accuracy will 
be severely degraded when integrating low-mass planets. 

Unfortunately, choosing a partition function that depends on velocity 
rather than position does not overcome this problem. A better way to tackle 
this issue is to use a set of coordinates that doesn't give rise to momentum 
cross terms like those in Hj ump or i?Small> so that indirect terms are a func- 
tion of the coordinates only. Barycentric coordinates have this property, but 
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as noted earlier, symplectic integrators that use barycentric coordinates tend 
to perform poorly. Jacobi coordinates won't work since they do not have the 
necessary properties for treating close encounters between planets. 

Chambers (2003) described a new set of coordinates with the right prop- 
erties to apply Levison and Duncan's scheme for integrating close encounters 
with a star. These coordinates, dubbed "Yosemite coordinates" by the author, 
are given by 



X 



/ moxo + Ej m i x j 

V m tot 



Xj = (xj - x ) + — ^2 rrijixij - x ) (52) 
mo . 

for the case of planets orbiting a single star, where the subscript refers to 
the star, and 

where 



H = 3 (54) 



m 

The canonically conjugate momenta are 

N 

Po = Po + Pj 

Po + (l + /3 + /3 M )E J P J 



m to t 



(55) 



Using these coordinates, the Hamiltonian for a system of planets orbiting 
a single star is 

H = H Kcp + H lBt (56) 

where 

A / if Gm m, 



i=l 



2m, i?. 



where 



_ ^EjmjXj 



m (l + /?M) (58) 
The second set of terms in iii nt represent indirect perturbations on each 
planet. These terms can be divided between iii nt and i?Kep in an analogous 
manner to the scheme of Levison et al. (2000) as follows: 
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_ / Pi _ Gm m t 

^Large ~ 2^\n ms 



i = l 
N 



2mi Ri 



#Small = - 2^ ' 



i— 1 J>2 ^ 



- / 1 1 \ 

£ G m „m, — - — — [1 - A(R U R 2 , R N )} (59) 



Ri X, — S 

2=1 

where A is chosen so that A — ► 1 whenever any planet approaches the star 
and j4 = when all planets are far from the star. 

The advantage of using Yosemite coordinates now becomes clear: each 
of the partitioned indirect terms in (|59|) is proportional to mj. As a result, 
low-mass planets and test particles can be integrated to the same accuracy as 
massive planets without the problems encountered when using (|49[) . The same 
scheme can be extended to binary systems by using new coordinate systems 
analogous to Yosemite coordinates with the binary hierarchy built in, as in 
the wide-binary and close-binary coordinate systems. 



6 Summary 

Conventional integration algorithms such as Runge-Kutta and Bulirsch-Stoer 
can be applied to planetary systems orbiting binary stars with little or no 
modification. However, these algorithms tend to be slow and exhibit long- 
term growth in energy errors. For these reasons, symplectic integration algo- 
rithms have become the tool of choice for many researchers. In this chapter, 
I have described how the standard mixed-variable symplectic map (MVS) of 
Wisdom and Holman (1991) can be adapted for use with planets in binary 
systems, including modifications to handle close approaches between planets 
with each other and with the stars themselves. In addition, I have shown that 
the performance of these algorithms can be improved at little extra cost by 
using symplectic correctors. These adaptations mean that symplectic algo- 
rithms can now be applied to a wide variety of problems involving planets in 
binary-star systems. However, there is still scope for future improvements. In 
particular, the current method for following close encounters between a planet 
and a star is slow and cumbersome, and there remains no easy way to handle 
planets whose orbital motion switches from one star to the other. 
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